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Abstract 

The variational determination of the two-particle density matrix is an interesting, but not yet 
fully explored technique that allows to obtain ground-state properties of a quantum many-body 
system without reference to an A^-particle wave function. The one-dimensional fermionic Hubbard 
model has been studied before with this method, using standard two- and three-index conditions 
on the density matrix [J. R. Hammond et al, Phys. Rev. A 73, 062505 (2006)], while a more 
recent study explored so-called subsystem constraints [N. Shenvi et al, Phys. Rev. Lett. 105, 
213003 (2010)]. These studies reported good results even with only standard two-index conditions, 
but have always been limited to the half- filled lattice. In this Letter we establish the fact that 
the two-index approach fails for other fillings. In this case, a subset of three-index conditions 
is absolutely needed to describe the correct physics in the strong- repulsion limit. We show that 
applying lifting conditions [J.R. Hammond et al., Phys. Rev. A 71, 062503 (2005)] is the most 
economical way to achieve this, while still avoiding the computationally much heavier three- index 
conditions. A further extension to spin-adapted lifting conditions leads to increased accuracy in 
the intermediate repulsion regime. At the same time we establish the feasibility of such studies to 
the more complicated phase diagram in two-dimensional Hubbard models. 

PACS numbers: 
Keywords: 



1 



The main problem in many-body quantum mechanics, which comprises nuclear physics, 
quantum chemistry and condensed matter physics, is the exponential increase of the di- 
mension of Hilbert space with the number of particles. The challenge has therefore been 
to develop approximate methods which describe the relevant degrees of freedom in the sys- 
tem without an excessive computational cost, i.e. with a polynomial increase. In one of 
these methods, the A^-particle wave function is replaced by the two-particle density matrix 
(2DM), and over the last decade, a lot of progress has been made in this field [1]-|6]. For a 
Hamiltonian: 

H = ^ ta^aicifi + ^ XI ^al3;isO'Wf}a'Sa'r > (1) 
containing only pairwise interactions, the energy of the system can be expressed as: 

EiV) = Tr ri7(2) = 1 Yl ^^P-nsH^'L ' (2) 

a/375 

in terms of the 2DM: 

Tap-ns = {^''\aiala5a^\m'') , (3) 
and the reduced two-particle Hamiltonian, 

(2) 

-^a/3;75 ~ N — 1 ^^""^^i^^ ~ ^aSip-f — SfSjtaS + Spsta-y) + VaPiyS ■ (4) 

Second-quantized notation is used where aj^ (aa) creates (annihilates) a fermion in the single- 
particle state a. 

In variational density-matrix optimization (v2DM), originally introduced by Lowdin, 
Mayer and Coleman [3^3, one exploits this fact and uses the 2DM as a variable in a vari- 
ational approach. From the resulting 2DM all one- and two-body properties of the ground 
state can be extracted. This should not be implemented naively, however, as there are a 
number of non-trivial constraints which a 2DM has to fulfil in order to be derivable from a 
A^-particle wave function. This is the A^-representability problem [9] which was proven to 
belong, in general, to the QMA-complete complexity class [10\. In practical approaches one 
uses a set of conditions which are necessary but not sufficient, and therefore lead to a lower 
bound on the ground-state energy. The most commonly used are the two-index conditions, 
called P (or D), Q and G [HI HI], and the computationally much heavier three-index con- 
ditions called Ti and T2 [121 113] • They all rely on the fact that for a manifestly positive 
Hamiltonian H = Yli^l^iy expectation value of the energy has to be larger than zero. 
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u 



PQG PQGT V2.5DM exact 



50 



-3.55 -2.29 



-2.28 -2.20 



100 



-3.49 -2.15 



-2.14 -2.08 



1000 



-3.44 -2.03 



-2.01 -2.01 



TABLE I: Ground-state energy of a 6-site lattice with 5 particles for U = 50, 100 and 1000, exact 
results compared with v2DM results using PQG and PQGT results, and v2.5DM results. 

These conditions can be expressed as linear matrix maps of the 2DM that have to be positive 
semidefinite. Another type of constraint that has recently been developed are the subsystem 
or active-space constraints [HHIS] in which linear conditions are imposed only that part of 
the density matrix that is related to a subspace of the complete single-particle space. This 
allows to increase accuracy (in the subspace) without having to use three-index conditions. 
Such v2DM methods have been used to study a wide variety of many-body systems: nuclei 
[IT] , atoms and molecules P-E], but also lattice systems [H [lH HEl [18], [19]. 

The Hubbard Hamiltonian [2Qj is the simplest schematic Hamiltonian that models the 
non-trivial correlations in solids as a competition between a delocalizing hopping term and 
an on-site repulsion term. In one dimension this Hamiltonian reads: 



where the sites on a periodic lattice are labeled a and a is the (up or down) spin. In 
previous v2DM studies of the one-dimensional Hubbard model [H [161 HH] only the half-filled 
lattice was studied, and it was found that even the two-index conditions could accurately 
describe the ground-state energy. In this Letter we show that the two-index conditions fail 
to describe the strong correlation limit below half filling, and that the subsystem constraints, 
as introduced in 116] . cannot solve this problem. In fact, a particular type of three-index 
conditions are needed in this limit, and we show that by using the 2.5DM (which is the 3DM 
diagonal in one spatial index) as the central object, these constraints can be incorporated 
while keeping the basic matrix manipulations in two-particle space. 

In order to demonstrate the problem. Fig. [T] shows the ground-state energy as a function 
of the on-site repulsion f/, for 5 particles in a 6-site lattice. As one can see the general form 
of the exact E vs. U curve (obtained through diagonalization) is nicely described by the 
PQGT result {i.e. two- and three-index conditions). On the other hand the PQG result 





3 




FIG. 1: Ground-state energy as a function of on-site repulsion U for 5 particles on a 6-site lattice. 
Exact results compared with v2DM results using PQG and PQGT1T2, and the v2.5DM results. 

(only two-index conditions) grossly underestimates the energy when U increases. The large- 
U limit is examined in Table [T] and one notes that PQG fails to get the energy right in this 
limit, as opposed to the PQGT result. When inspecting the PQG-optimized 2DM in the 
large-?/ limit, it was found that the on-site repulsion term vanishes, as it should be, since 
the 2DM elements corresponding to doubly occupied sites are zero: 

So the problem with the two-index conditions lies in its inability to describe the hopping 
term on a lattice where the sites cannot be doubly occupied. It is readily understood that 
subsystem constraints cannot solve this issue, because the singly-occupied space is a subspace 
of the full A^-particle Hilbert space, and cannot be obtained by a restriction of single-particle 
space to a subsystem. 

The creation and annihilation of particles on a singly occupied lattice can be described 
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by the so-called Gutzwiller operators [211 [22] : 



t \ t 

1 CL^OiQi j Qj^ , 



(7) 
(8) 



where a and a are single particle indices on the same site with opposite spin. In analogy 
with the necessary and sufficient conditions for A^-representability of the one-particle density 
matrix (IDM) [H], one can state that all Hamiltonians expressed as ffist-order operators of 
the (7a's will be correctly optimized if the following 'Gutzwiller' matrix positivity conditions 
are satisfied: 



^ 



with p% = {^''\gig,\^'') 

^ with q% = {^''\gj^\^'') 



(9) 
(10) 



These matrices can be expressed as a function of regular fermionic creation and annihilation 

operators, 

P% = Pa(S - ^a^-M ~ ^aa;l3a + (^^ | alaUa^a/fa/J | , (H) 

(12) 



QaP — ^a/B (1 — Paa — P/J/?) + ^ I3l3;l3a + ^al3;aa + | Ct^aaaQa^Cl lO/? | ) 



in which p is the IDM: 



p^p = {^'^\aiap\^'') 



(13) 



It is clear from Eqs. (11), (12) that the 3DM plays an essential role in describing the strong 



correlation limit (which is why the two-index conditions fail) but one also sees that not the 
full three-body space is needed. In fact, the 2.5DM, defined as: 



W 



l\SiSab]Scd) 

\ab;cd 



2 C _|_ ]^ / ^ ' I \abl Icdl I / 

M 



(14) 



is the minimal object from which both the PQG conditions and the Gutzwiller conditions 
can be derived, and for which basic matrix manipulations are still on two-particle space. In 



Eq. (14) creates three particles: 

-° labl 



Sab 



® a 



(15) 



on lattice sites a, b and /, coupled to total spin S, spin projection Sz and intermediary spin 
Sab- Note that the spin-averaged ensemble is used [5] for describing the A^-particle state 



with total spin S. Here one considers an equal weight ensemble of all spin projections A^, 
and as a result the 2.5DM has no Sz dependence. The 2.5DM is a block-diagonal object, 
in the sense that it is the 3DM diagonal in one pair of spatial indices. It can be used as 
the central object in a variational approach, applying constraints that include the PQG and 
Gutzwiller conditions. This is a generalization of an approach used by Mazziotti et al. [TlfTS] 
in a discussion of partial three-positivity constraints. They used a 3DM diagonal in both 
spatial and spin indices as a variational object in a study of the Lipkin spin model. Letting 
the spin index be off-diagonal allows us to construct a spin-coupled version of the 2.5DM, 
which leads to an increase in speed of the optimization (see e.g. jS]). Apart from this, the 
increased flexibility of the 2.5DM is important as it captures more correlation, which leads 
to a better result for the ground-state energy. We find that going from the spin-uncoupled 
to the spin-coupled form removes about 20% from the remaining discrepancy with the exact 
result, in the intermediate U/t region 

The first non-trivial constraint we impose on the 2.5DM is a consistency condition that 
ensures symmetry between the diagonal third index and the other indices. As an example, 
one of these relations reads: 

W'^'"> = |S„J|SJ Y: \S.A\SA If i ''"I If \ '"]w>^''"> , (16) 

with [5*] = \/2S + 1. In addition, we add constraints that are analogous to the standard 
two- and three-index conditions in that they can be expressed as matrix maps of the 2.5DM 
that have to be positive semidefinite. The first condition is simply that the different blocks 
of W have to be positive semidefinite: 

ly' ^ . (17) 

The other five conditions are spin-adapted generalizations of the lifting conditions introduced 
in [H [131 E3] , and are of the form: 

C {W)' >z with C {Wf iS^^-) = E(^^^l^^l2f"^^l2f"^l^^^) , (18) 

M 

in which the -B^ consist of different combinations of creation and annihilation operators. As 
an example of such a condition, consider 5^ defined as: 

s \ 

) where a^^^ = (-l)^+"''^aa_^^ . (19) 

MabJ 
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FIG. 2: Ground-state energy as a function of on-site repulsion U for 9 particles on a 10-site lattice. 
Exact results compared with v2.5DM results using 2.5-index conditions and v2DM results using 
PQG. 

The various £'s arise by considering (aaa), (a^a^a), (aaa"'^), [aa^a^) and (a^aa) combinations, 
and can all be expressed as a function of W through the use of anticommutation relations 
and spin recouphng. 

The numerical optimization of the 2.5DM under these positivity constraints is a semidef- 
inite program, and exactly the same methods used for the optimization of the 2DM can 
be used [H [211 125] • The scaling of the basic matrix manipulations in this optimization is 
M"^, as opposed to the full three-index conditions, which scale as M^, with M the size of 
single-particle Hilbert space. The result of such a v2.5DM calculation is also shown in Fig. [T} 
It is clear that, as anticipated, the strong interaction limit is described with PQGT quality 
without resorting to the full PQGT framework. In fact, the v2.5DM results are slightly 
better then the PQGT. This is because the Ti and T2 conditions express the positivity 
of an anticommutator of three-particle operators, whereas in v2.5DM positivity is imposed 
on all possible individual products of three-particle operators, be it of a restricted class. 
Similar results are obtained for a somewhat larger lattice of 10 sites, shown in Fig. [2j In this 
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u 



PQG PQGT V2.5DM exact 



50 



-4.89 -2.54 



-2.53 -2.46 



100 



-4.77 -2.27 



-2.26 -2.22 



1000 



-4.67 -2.03 



-2.03 -2.02 



TABLE II: Ground-state energy of a 10-site lattice with 9 particles for U = 50, 100 and 1000, 
exact results compared with v2DM results using PQG and PQGT results, and v2.5DM results. 

case direct diagonalization is no longer an option, but we can compare to quasi-exact results 
calculated with a matrix product state (MPS) optimization [26H28|, which is uniquely suited 
for this kind of one-dimensional lattice problem. Again, the v2DM result using PQG con- 
ditions is inaccurate, and one has to incorporate the three-particle correlations captured in 
the 2.5DM approach. The full-blown PQGT calculation is far more costly than v2.5DM but 
produces slightly inferior results. Clearly both the [/ — )■ and strongly interacting U ^ oo 
limit are now exact. The latter statement is demonstrated in Tables U and [TTl and follows 
from the above discussion about the Gutzwiller conditions. 

In summary, we developed the v2.5DM method that takes into account the necessary 
correlations needed to describe the large-f/ limit of the Hubbard model, without having to 
resort to full-blown three-index conditions. It must be stressed that up to know we have 
only included the spin symmetry of the model in our code. If translational symmetry, parity 
and pseudospin symmetry are taken into account much larger lattices can be considered. 
As an example, our fully symmetric PQG version allows lattice sizes up to 100 sites, and 
the fully symmetric PQGT program up to 20 sites. We expect a fully symmetric version 
of v2.5DM to be applicable to lattice sizes of about 50 sites, thereby enabling us to study 
two-dimensional lattices of reasonable size. 

The diagonality of the third index in the 2.5DM implies that the result will depend on 
the chosen single-particle basis. For the Hubbard model it is clear that the site basis is 
the optimal basis to use for the diagonal third index. It will be interesting to study other 
systems where it is less clear what the best choice of the single-particle basis would be. An 
appealing application, e.g., are molecules, where one can hope to get three-index (T1T2) 
precision by applying v2.5DM method with a carefully chosen basis. A first guess of what 
the best basis would be is the basis of natural orbitals, for which it has been shown that the 



8 



fuU-CI expansion has the fastest convergence [7] . This should be implemented using a outer 
self-consistency loop as a natural basis changes after each 2.5DM variational calculation. 
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